Exploratory Analysis for MoTR Reading Data

shhh <- suppressPackageStartupMessages # It's a library, so shhh!

shhh(library( mgcv ))
shhh(library(dplyr))
shhh(library(ggplot2))
shhh(library(lme4))
shhh(library(tidymv))
shhh(library(gamlss))
shhh(library(gsubfn))

Warning: unable to load shared object '/Library/Frameworks/R.framework/Resources/modules//R_X11.so':
  dlopen(/Library/Frameworks/R.framework/Resources/modules//R_X11.so, 0x0006): Library not loaded: '/opt/X11/lib/libSM.6.dylib'
  Referenced from: '/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/modules/R_X11.so'
  Reason: tried: '/opt/X11/lib/libSM.6.dylib' (no such file), '/usr/local/lib/libSM.6.dylib' (no such file), '/usr/lib/libSM.6.dylib' (no such file)tcltk DLL is linked to '/opt/X11/lib/libX11.6.dylib'

shhh(library(lmerTest))
shhh(library(tidyverse))
shhh(library(boot))
shhh(library(rsample))
shhh(library(plotrix))
shhh(library(ggrepel))
shhh(library(mgcv))

theme_set(theme_bw())
options(digits=4)
options(scipen=999)
set.seed(444)
pipe_message = function(.data, status) {message(status); .data}

Read in MoTR Data

rate = 160

file_prefix = "/Users/cui/Desktop/motr_provo_160/"
fnames = list.files(path=file_prefix)

df = data.frame()
for (f in fnames) {
  temp = read.csv(paste0(file_prefix, "/", f)) %>%
    mutate(subj = str_remove(f, "_reading_measures.csv")) %>%
    rename(go_past_time = go_pass_time)
  df = rbind(df, temp)
}

filter_df = df %>%
  group_by(para_nr, subj) %>%
    summarise(correct = if_else(unique(correctness) == 1, 1, 0)) %>%
  ungroup() %>%
  drop_na() %>%
  group_by(subj) %>%
    summarise(p_correct = mean(correct)) %>%
  ungroup() %>%
  mutate(p_correct = round(p_correct, digits = 2))

`summarise()` has grouped output by 'para_nr'. You can override using the `.groups` argument.

filter_df = filter_df %>%
  filter(p_correct < 0.8)
# View(filter_df)
## reader_3:0.70, reader_60:0.79, reader_76:0.72 , reader_256:0.71 , reader_262:0.57 

raw_df = df %>%
  mutate(word = str_trim(word)) %>%
  mutate(subj = str_remove(subj, "reader_")) %>%
  mutate(subj = as.integer(subj)) %>%
  filter(! subj %in% c(3, 60, 76, 256, 262)) %>%
  # See below for filtering out reading measures that are super long
  dplyr::select(expr_id, cond_id, para_nr, word, word_nr, first_duration, total_duration, gaze_duration, go_past_time, FPReg, subj) %>%
  drop_na()

raw_df

# Average across subjects
motr_agg_df = raw_df %>%
  gather(metric, value, 6:10) %>%
    group_by(para_nr, word_nr, word, metric) %>%
    # for non-binary measures, filter out that are super long
    mutate(outlier = if_else(metric != "FPReg" & value > (mean(value) + 3 * sd(value) ), T, F)) %>%
    filter(outlier == F) %>%
  # # Filter out words with a reading-time of zero
  # mutate(zero = if_else(metric != "FPReg" & value == 0, T, F)) %>%
  # filter(zero == F) %>%
  drop_na() %>%
    summarise(value = mean(value),
              nsubj = length(unique(subj))) %>%
  ungroup() %>%
  arrange(para_nr, word_nr) %>%
  rename(
    text_id = para_nr,
    word_text_idx = word_nr,
    motr_value = value
  )

`summarise()` has grouped output by 'para_nr', 'word_nr', 'word'. You can override using the `.groups` argument.

motr_agg_df

NA

Comparison to Provo

# Read in Provo surprisal, frequency and length data
provo_modeling_df = read.csv("/Users/cui/Documents/ETH/MoTR/pipeline/ancillary_data/provo_df.csv") %>%
  dplyr::select(text_id, sent_id, trigger_idx, word, freq, surp, len) %>%
  rename(word_idx = trigger_idx)

provo_modeling_df

NA

# Read in Provo eyetracking data

provo_raw_df = read.csv("/Users/cui/Documents/ETH/MoTR/pipeline/ancillary_data/provo_eyetracking.csv")

provo_eyetracking_df = provo_raw_df %>%
  dplyr::select(Participant_ID, Text_ID, Sentence_Number, Word_In_Sentence_Number, Word, Word_Number, IA_FIRST_FIX_PROGRESSIVE, IA_FIRST_RUN_DWELL_TIME, IA_DWELL_TIME, IA_REGRESSION_PATH_DURATION, IA_REGRESSION_OUT) %>%
  rename( #first_duration = IA_FIRST_FIXATION_DURATION,   
          gaze_duration = IA_FIRST_RUN_DWELL_TIME,
          total_duration = IA_DWELL_TIME,
          go_past_time = IA_REGRESSION_PATH_DURATION,
          subj = Participant_ID,
          text_id = Text_ID,
          sent_id = Sentence_Number,
          word_idx = Word_In_Sentence_Number,
          word_text_idx = Word_Number,   # IA_ID?
          word = Word,      
          FPReg = IA_REGRESSION_OUT,
          ff_progressive = IA_FIRST_FIX_PROGRESSIVE) %>%
  mutate(first_duration = gaze_duration) %>%
  mutate(gaze_duration = if_else(ff_progressive == 0, 0, gaze_duration),
         go_past_time = if_else(ff_progressive == 0, 0, go_past_time)) %>%
  dplyr::select(-ff_progressive) %>%
  gather(metric, value, 7:11) %>%
  mutate(value = if_else(is.na(value), as.integer(0), as.integer(value))) %>%
  drop_na() %>%
  mutate(word = str_trim(word)) %>%
  mutate(subj = str_remove(subj, "Sub")) %>%
  mutate(subj = as.integer(subj)) %>%
    group_by(text_id, word_text_idx, sent_id, word_idx, word, metric) %>%
    mutate(outlier = if_else(metric != "FPReg" & value > (mean(value) + 3 * sd(value) ), T, F)) %>%
    filter(outlier == F) %>%
  ungroup() #%>%
  # # Filter out words with a reading-time of zero
  # mutate(zero = if_else(metric != "FPReg" & value == 0,T, F)) %>%
  # filter(zero == F)

# Aggregate cross-participant data for all subjects
provo_eyetracking_agg_df = provo_eyetracking_df %>%
  group_by(text_id, word_text_idx, sent_id, word_idx, word, metric) %>%
    summarise(value = mean(value), .groups = 'drop') 

# Split the eyetracking data in two by subjects to see how well it correlates with itself
provo_eyetracking_subj1_df = provo_eyetracking_df %>%
  filter(subj <= 42) %>%
  group_by(text_id, word_text_idx, sent_id, word_idx, word, metric) %>%
    summarise(value = mean(value)) %>%
  ungroup() %>%
  rename(value_1 = value) %>%
  dplyr::select(-sent_id, -word_idx)

`summarise()` has grouped output by 'text_id', 'word_text_idx', 'sent_id', 'word_idx', 'word'. You can override using the `.groups` argument.

  
  provo_eyetracking_subj2_df = provo_eyetracking_df %>%
  filter(subj > 42) %>%
  group_by(text_id, word_text_idx, sent_id, word_idx, word, metric) %>%
    summarise(value = mean(value)) %>%
  ungroup() %>%
    rename(value_2 = value)%>%
  dplyr::select(-sent_id, -word_idx)

`summarise()` has grouped output by 'text_id', 'word_text_idx', 'sent_id', 'word_idx', 'word'. You can override using the `.groups` argument.

  
provo_eyetr_grouped_df = merge(provo_eyetracking_subj2_df, provo_eyetracking_subj1_df, by=c("text_id", "word_text_idx", "metric")) %>%
  filter(word.x == word.y) %>%
  dplyr::select(-word.x, -word.y) %>%
  group_by(metric) %>%
    mutate(motr_outlier = if_else(metric != "FPReg" & value_1 > (mean(value_1) + 3 * sd(value_1) ), T, F)) %>%
    filter(motr_outlier == F) %>%
    mutate(eyetr_outlier = if_else(metric != "FPReg" & value_2 > (mean(value_2) + 3 * sd(value_2) ), T, F)) %>%
    filter(eyetr_outlier == F) %>%
  ungroup() %>%
  gather(measure, value, c("value_1", "value_2")) %>%
  dplyr::select(-motr_outlier, -eyetr_outlier)

provo_df = merge(provo_eyetracking_agg_df, provo_modeling_df, by=c("text_id", "sent_id", "word_idx")) %>%
  mutate(word_text_idx = as.integer(word_text_idx - 1)) %>%
  arrange(text_id, sent_id, word_idx) %>%
  rename(eyetr_value = value) 

provo_df = merge(provo_df, motr_agg_df, by=c("text_id", "word_text_idx", "metric")) %>%
arrange(text_id, sent_id, word_idx) %>%
  # almost all the word.x != word.y is because of normalization problem, so we can keep them, instead, deleting some special cases
  filter(!(text_id == 13 & word_text_idx >= 20 & word_text_idx <= 52)) %>%
  filter(!(text_id == 3 & word_text_idx >= 46 & word_text_idx <= 57)) %>%
# # filter(word.y == word) %>%
dplyr::select(-word.x, -word.y) %>%
group_by(metric) %>%
  mutate(motr_outlier = if_else(metric != "FPReg" & motr_value > (mean(motr_value) + 3 * sd(motr_value) ), T, F)) %>%
  filter(motr_outlier == F) %>%
  mutate(eyetr_outlier = if_else(metric != "FPReg" & eyetr_value > (mean(eyetr_value) + 3 * sd(eyetr_value) ), T, F)) %>%
  filter(eyetr_outlier == F) %>%
ungroup() %>%
gather(measure, value, c("eyetr_value", "motr_value")) %>%
  dplyr::select(-motr_outlier, -eyetr_outlier)

provo_df %>%
  mutate(measure = if_else(measure == "eyetr_value", "Eyetracking Value", "MoTR Value")) %>%
  filter(metric != "FPReg") %>%
  ggplot(aes(x = value, color=metric)) +
    geom_density() +
    facet_wrap(.~measure, scales="free_y") +
    xlab("Reading Time (ms)")

# ggsave("../visualization/density.png", device = "png", width = 6, height = 2.5)

print("Gaze Duration")

[1] "Gaze Duration"

cor_df = provo_df %>% filter(metric == "gaze_duration") %>% spread(measure, value)
print(cor.test(cor_df$eyetr_value, cor_df$motr_value)$estimate)

   cor 
0.7874 

cor_df = provo_eyetr_grouped_df %>% filter(metric == "gaze_duration") %>% spread(measure, value)
print(cor.test(cor_df$value_1, cor_df$value_2)$estimate)

   cor 
0.9168 

print("First Duration")

[1] "First Duration"

gd_df = provo_df %>% filter(metric == "first_duration") %>% spread(measure, value)
print(cor.test(gd_df$eyetr_value, gd_df$motr_value)$estimate)

   cor 
0.7815 

cor_df = provo_eyetr_grouped_df %>% filter(metric == "first_duration") %>% spread(measure, value)
print(cor.test(cor_df$value_1, cor_df$value_2)$estimate)

   cor 
0.9201 

print("Go Past Time")

[1] "Go Past Time"

gd_df = provo_df %>% filter(metric == "go_past_time") %>% spread(measure, value)
print(cor.test(gd_df$eyetr_value, gd_df$motr_value)$estimate)

   cor 
0.7292 

cor_df = provo_eyetr_grouped_df %>% filter(metric == "go_past_time") %>% spread(measure, value)
print(cor.test(cor_df$value_1, cor_df$value_2)$estimate)

   cor 
0.9133 

print("Total Duration")

[1] "Total Duration"

gd_df = provo_df %>% filter(metric == "total_duration") %>% spread(measure, value)
print(cor.test(gd_df$eyetr_value, gd_df$motr_value)$estimate)

   cor 
0.7601 

cor_df = provo_eyetr_grouped_df %>% filter(metric == "total_duration") %>% spread(measure, value)
print(cor.test(cor_df$value_1, cor_df$value_2)$estimate)

   cor 
0.9275 

print("Regression")

[1] "Regression"

reg_df = provo_df %>% filter(metric == "FPReg") %>% spread(measure, value)
print(cor.test(reg_df$eyetr_value, reg_df$motr_value)$estimate)

   cor 
0.3259 

cor_df = provo_eyetr_grouped_df %>% filter(metric == "FPReg") %>% group_by(text_id, metric, measure) %>%
  summarize(value = mean(value, na.rm = TRUE), .groups = 'drop') %>% spread(measure, value)
print(cor.test(cor_df$value_1, cor_df$value_2)$estimate)

   cor 
0.6807 

provo_df %>%
  filter(metric != "FPReg") %>%
  spread(measure, value) %>%
  ggplot(aes(x = motr_value, y=eyetr_value)) +
    geom_point(alpha = 0.05) +
    geom_abline(slope=1, intercept=0, color = "black") +
    #stat_summary_bin(bins=100, fun.data = "mean_cl_boot", size = 0.05) +
    facet_wrap(.~metric, scales = "free", nrow = 1) +
    coord_cartesian(ylim=c(0, 500), xlim=c(0, 500)) +
    geom_smooth()

# ggsave("../visualization/metric_cor.png", device = "png", width = 6, height = 2.5)

Correlations to Word-Level Statistical Properties

print("Gaze Duration")

[1] "Gaze Duration"

print("Len")

[1] "Len"

cor_df = provo_df %>% filter(metric == "gaze_duration") %>% spread(measure, value)
print(cor.test(cor_df$eyetr_value, cor_df$len)$estimate)

   cor 
0.8597 

print(cor.test(cor_df$motr_value, cor_df$len)$estimate)

   cor 
0.8644 

print("Freq")

[1] "Freq"

cor_df = provo_df %>% filter(metric == "gaze_duration") %>% spread(measure, value)
print(cor.test(cor_df$eyetr_value, cor_df$freq)$estimate)

    cor 
-0.8069 

print(cor.test(cor_df$motr_value, cor_df$freq)$estimate)

    cor 
-0.7454 

print("Surp")

[1] "Surp"

cor_df = provo_df %>% filter(metric == "gaze_duration") %>% spread(measure, value)
print(cor.test(cor_df$eyetr_value, cor_df$surp)$estimate)

   cor 
0.5683 

print(cor.test(cor_df$motr_value, cor_df$surp)$estimate)

   cor 
0.4978 

print("Total Duration")

[1] "Total Duration"

print("Len")

[1] "Len"

cor_df = provo_df %>% filter(metric == "total_duration") %>% spread(measure, value)
print(cor.test(cor_df$eyetr_value, cor_df$len)$estimate)

  cor 
0.825 

print(cor.test(cor_df$motr_value, cor_df$len)$estimate)

  cor 
0.838 

print("Freq")

[1] "Freq"

cor_df = provo_df %>% filter(metric == "total_duration") %>% spread(measure, value)
print(cor.test(cor_df$eyetr_value, cor_df$freq)$estimate)

   cor 
-0.781 

print(cor.test(cor_df$motr_value, cor_df$freq)$estimate)

    cor 
-0.7258 

print("Surp")

[1] "Surp"

cor_df = provo_df %>% filter(metric == "total_duration") %>% spread(measure, value)
print(cor.test(cor_df$eyetr_value, cor_df$surp)$estimate)

   cor 
0.5931 

print(cor.test(cor_df$motr_value, cor_df$surp)$estimate)

   cor 
0.5104 

print("First Duration")

[1] "First Duration"

print("Len")

[1] "Len"

cor_df = provo_df %>% filter(metric == "first_duration") %>% spread(measure, value)
print(cor.test(cor_df$eyetr_value, cor_df$len)$estimate)

   cor 
0.8444 

print(cor.test(cor_df$motr_value, cor_df$len)$estimate)

   cor 
0.8598 

print("Freq")

[1] "Freq"

cor_df = provo_df %>% filter(metric == "first_duration") %>% spread(measure, value)
print(cor.test(cor_df$eyetr_value, cor_df$freq)$estimate)

    cor 
-0.8016 

print(cor.test(cor_df$motr_value, cor_df$freq)$estimate)

    cor 
-0.7453 

print("Surp")

[1] "Surp"

cor_df = provo_df %>% filter(metric == "first_duration") %>% spread(measure, value)
print(cor.test(cor_df$eyetr_value, cor_df$surp)$estimate)

   cor 
0.5852 

print(cor.test(cor_df$motr_value, cor_df$surp)$estimate)

  cor 
0.507 

print("Go Past Time")

[1] "Go Past Time"

print("Len")

[1] "Len"

cor_df = provo_df %>% filter(metric == "go_past_time") %>% spread(measure, value)
print(cor.test(cor_df$eyetr_value, cor_df$len)$estimate)

   cor 
0.8063 

print(cor.test(cor_df$motr_value, cor_df$len)$estimate)

   cor 
0.8124 

print("Freq")

[1] "Freq"

cor_df = provo_df %>% filter(metric == "go_past_time") %>% spread(measure, value)
print(cor.test(cor_df$eyetr_value, cor_df$freq)$estimate)

    cor 
-0.7516 

print(cor.test(cor_df$motr_value, cor_df$freq)$estimate)

   cor 
-0.702 

print("Surp")

[1] "Surp"

cor_df = provo_df %>% filter(metric == "go_past_time") %>% spread(measure, value)
print(cor.test(cor_df$eyetr_value, cor_df$surp)$estimate)

   cor 
0.5476 

print(cor.test(cor_df$motr_value, cor_df$surp)$estimate)

   cor 
0.4893 

provo_df %>%
  gather(word_prop, word_prop_val, c("freq", "len", "surp")) %>%
  mutate(measure = if_else(measure == "eyetr_value", "Eyetracking Value", "MoTR Value")) %>%
  mutate(word_prop = case_when(
    word_prop == "freq" ~ "Frequency",
    word_prop == "len" ~ "Length",
    word_prop == "surp" ~ "Surprisal"
  )) %>%
  filter(metric == "gaze_duration") %>%
  ggplot(aes(x = value, y=word_prop_val, color = measure)) +
    geom_point(alpha = 0.1) +
    facet_wrap(measure~word_prop, scales="free", strip.position = "right") +
    geom_smooth(color = "grey") +
    xlab("Reading Measure") +
  theme(
    legend.position = "none",
    strip.placement = "outside"
  )

# ggsave("../visualization/word_prop_comps.png", device = "png", width = 6, height = 3)

provo_df %>%
  ggplot(aes(x = value, y=freq, color=metric)) +
    geom_point(alpha = 0.1) +
    facet_grid(metric~measure, scales="free") +
    geom_smooth()

provo_df %>%
  ggplot(aes(x = value, y=surp, color=metric)) +
    geom_point(alpha = 0.2) +
    facet_grid(metric~measure, scales="free") +
    geom_smooth()

Shape of surprisal / RT relationship

for current word:

fit_gam_inner = function(bootstrap_sample, mean_predictors) {
  
  df = bootstrap_sample$data
  weights = tabulate(as.integer(bootstrap_sample), nrow(df))
  
  m = gam(psychometric ~ s(surp, bs = 'cr', k = 6) + s(prev_surp, bs = 'cr', k = 6) + te(freq, len, bs = 'cr') + te(prev_freq, prev_len, bs = 'cr'), data = df, weights = weights)
  terms_to_predict = c("s(surp)", "s(prev_surp)")
  
  newdata = data.frame(surp=seq(0,20,by=0.1), prev_surp=mean_predictors$surp,
                       #surp=mean_predictors$surp, prev_surp=seq(0,20,by=0.1),
                       freq=mean_predictors$freq, prev_freq=mean_predictors$freq,
                       len=mean_predictors$len, prev_len=mean_predictors$len)

  # Returns a matrix N_samples * N_terms.
  per_term_predictions = predict(m, newdata=newdata, terms=terms_to_predict, type="terms")

  # Additive model -- sum across predictor response contributions (matrix columns).
  predictions = rowSums(per_term_predictions)

  return(newdata %>% mutate(y=predictions))
}

fit_gam = function(df, mean_predictors, alpha=0.05) {
  # Bootstrap-resample data
  boot_models = df %>% bootstraps(times=10) %>% 
   # Fit a GAM and get predictions for each sample
    mutate(smoothed=map(splits, fit_gam_inner, mean_predictors=mean_predictors))
  
  # Extract mean and 5% and 95% percentile y-values for each surprisal value
  result = boot_models %>% 
    unnest(smoothed) %>% 
    dplyr::select(surp, y) %>% 
    #dplyr::select(prev_surp, y) %>% 
    group_by(surp) %>% 
    #group_by(prev_surp) %>%
      summarise(y_lower=quantile(y, alpha / 2), 
                y_upper=quantile(y, 1 - alpha / 2),
                y=mean(y)) %>% 
    ungroup()
  
  return (result)
}

gam_modeling_df = provo_df %>%
  spread(measure, value) %>%
  # mutate(len = nchar(word)) %>%      # len has already exists, but do not count punct into len.
  group_by(metric, text_id) %>%
    arrange(word_text_idx) %>%
    mutate(prev_surp = lag(surp),
           prev_freq = lag(freq),
           prev_len = lag(len),
           prev_eyetr_value = lag(eyetr_value)) %>%
  ungroup() %>%
  drop_na() %>%
  rename(psychometric = motr_value)

smooths_df = data.frame()

metrics = c("gaze_duration", "total_duration", "go_past_time", "first_duration")
for (m in metrics) {
  print(paste0("Fitting model for ", m))
  dummy_df = gam_modeling_df %>% filter(metric == m)
  mean_predictors = dummy_df %>% summarise(surp = mean(surp), len = mean(len), freq = mean(freq))
  smooths = dummy_df %>% fit_gam(., mean_predictors)
  #Fix 0 surprisal = 0 ms
  gam_smooths = smooths %>% mutate(delta = 0 - y[1], y=y + delta, y_lower= y_lower + delta, y_upper=y_upper + delta)
  smooths_df = rbind(smooths_df, gam_smooths %>% mutate(psychometric = m))
}

[1] "Fitting model for gaze_duration"
[1] "Fitting model for total_duration"
[1] "Fitting model for go_past_time"
[1] "Fitting model for first_duration"

Get Density Data

get_d_points = function(df) {
    x = density(df$surp)$x
    y = density(df$surp)$y
    return(data.frame(x, y))
  }

density_data = data.frame()

for(m in c("gaze_duration", "total_duration", "go_past_time", "first_duration")) {
  dummy_df = provo_df %>% filter(metric == m) %>%
      do({get_d_points(.)}) %>%
      filter(x>0, x<20)
  density_data = rbind(density_data, dummy_df %>% mutate(metric=m))
}

# Surprisal curves
  ggplot() +
      # Density Data
      annotate("rect", xmin=0, xmax=20, ymin=-20,ymax=-10, fill="#f4f4f4", color="grey", alpha=1, size = 0) +
      geom_line(data = density_data, aes(x=x, y=y*50 - 18), color="#aaaaaa", size = 0.4) +
      # Surrp / Rt data
      #geom_line(data = smooths_df, aes(x=prev_surp, y=y, color = psychometric), size=0.7) +
      geom_line(data = smooths_df, aes(x=surp, y=y, color = psychometric), size=0.7) +
      #geom_ribbon(data = smooths_df, aes(x=prev_surp, ymin=y_lower, ymax=y_upper, fill = psychometric), alpha=0.3, size=0.5) +
      geom_ribbon(data = smooths_df, aes(x=surp, ymin=y_lower, ymax=y_upper, fill = psychometric), alpha=0.3, size=0.5) +
      scale_x_continuous(labels=c(0, 10, 20), breaks=c(0, 10, 20), minor_breaks = NULL) +
      facet_wrap(psychometric~., nrow = 1) +
      ylab("Slowdown due to Surprisal (ms)") +
      xlab("Surprisal of Word") +
      ggtitle("MoTR Times and Current Word Surprisal")

  theme(
    legend.position = "none",
    panel.grid.minor = element_blank()
  )

List of 2
 $ legend.position : chr "none"
 $ panel.grid.minor: list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 - attr(*, "class")= chr [1:2] "theme" "gg"
 - attr(*, "complete")= logi FALSE
 - attr(*, "validate")= logi TRUE

  
# ggsave("../visualization/surprisal_rt_link.png", device = "png", width = 6, height = 2.5)

for previous word:

fit_gam_inner_2 = function(bootstrap_sample, mean_predictors) {
  
  df = bootstrap_sample$data
  weights = tabulate(as.integer(bootstrap_sample), nrow(df))
  
  m = gam(psychometric ~ s(surp, bs = 'cr', k = 6) + s(prev_surp, bs = 'cr', k = 6) + te(freq, len, bs = 'cr') + te(prev_freq, prev_len, bs = 'cr'), data = df, weights = weights)
  terms_to_predict = c("s(surp)", "s(prev_surp)")
  
  newdata = data.frame(surp=mean_predictors$surp, prev_surp=seq(0,20,by=0.1),
                       freq=mean_predictors$freq, prev_freq=mean_predictors$freq,
                       len=mean_predictors$len, prev_len=mean_predictors$len)

  # Returns a matrix N_samples * N_terms.
  per_term_predictions = predict(m, newdata=newdata, terms=terms_to_predict, type="terms")

  # Additive model -- sum across predictor response contributions (matrix columns).
  predictions = rowSums(per_term_predictions)

  return(newdata %>% mutate(y=predictions))
}

fit_gam_2 = function(df, mean_predictors, alpha=0.05) {
  # Bootstrap-resample data
  boot_models = df %>% bootstraps(times=10) %>% 
   # Fit a GAM and get predictions for each sample
    mutate(smoothed=map(splits, fit_gam_inner_2, mean_predictors=mean_predictors))
  
  # Extract mean and 5% and 95% percentile y-values for each surprisal value
  result = boot_models %>% 
    unnest(smoothed) %>% 
    dplyr::select(prev_surp, y) %>%
    group_by(prev_surp) %>%
      summarise(y_lower=quantile(y, alpha / 2), 
                y_upper=quantile(y, 1 - alpha / 2),
                y=mean(y)) %>% 
    ungroup()
  
  return (result)
}

gam_modeling_df_2 = provo_df %>%
  spread(measure, value) %>%
  # mutate(len = nchar(word)) %>%  # len has already exists, but do not count punct into len.
  group_by(metric, text_id) %>%
    arrange(word_text_idx) %>%
    mutate(prev_surp = lag(surp),
           prev_freq = lag(freq),
           prev_len = lag(len),
           prev_eyetr_value = lag(eyetr_value)) %>%
  ungroup() %>%
  drop_na() %>%
  rename(psychometric = motr_value)

smooths_df = data.frame()

metrics = c("gaze_duration", "total_duration", "go_past_time", "first_duration")
for (m in metrics) {
  print(paste0("Fitting model for ", m))
  dummy_df = gam_modeling_df_2 %>% filter(metric == m)
  mean_predictors = dummy_df %>% summarise(surp = mean(surp), len = mean(len), freq = mean(freq))
  smooths = dummy_df %>% fit_gam_2(., mean_predictors)
  #Fix 0 surprisal = 0 ms
  gam_smooths = smooths %>% mutate(delta = 0 - y[1], y=y + delta, y_lower= y_lower + delta, y_upper=y_upper + delta)
  smooths_df = rbind(smooths_df, gam_smooths %>% mutate(psychometric = m))
}

[1] "Fitting model for gaze_duration"
[1] "Fitting model for total_duration"
[1] "Fitting model for go_past_time"
[1] "Fitting model for first_duration"

get_d_points = function(df) {
    x = density(df$surp)$x
    y = density(df$surp)$y
    return(data.frame(x, y))
  }

density_data = data.frame()

for(m in c("gaze_duration", "total_duration", "go_past_time", "first_duration")) {
  dummy_df = provo_df %>% filter(metric == m) %>%
      do({get_d_points(.)}) %>%
      filter(x>0, x<20)
  density_data = rbind(density_data, dummy_df %>% mutate(metric=m))
}

# Surprisal curves
  ggplot() +
      # Density Data
      annotate("rect", xmin=0, xmax=20, ymin=-20,ymax=-10, fill="#f4f4f4", color="grey", alpha=1, size = 0) +
      geom_line(data = density_data, aes(x=x, y=y*50 - 18), color="#aaaaaa", size = 0.4) +
      # Surrp / Rt data
      geom_line(data = smooths_df, aes(x=prev_surp, y=y, color = psychometric), size=0.7) +
      geom_ribbon(data = smooths_df, aes(x=prev_surp, ymin=y_lower, ymax=y_upper, fill = psychometric), alpha=0.3, size=0.5) +
      scale_x_continuous(labels=c(0, 10, 20), breaks=c(0, 10, 20), minor_breaks = NULL) +
      facet_wrap(psychometric~., nrow = 1) +
      ylab("Slowdown due to Surprisal (ms)") +
      xlab("Surprisal of Word") +
      ggtitle("MoTR Times and Previous Word Surprisal")

  theme(
    legend.position = "none",
    panel.grid.minor = element_blank()
  )

List of 2
 $ legend.position : chr "none"
 $ panel.grid.minor: list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 - attr(*, "class")= chr [1:2] "theme" "gg"
 - attr(*, "complete")= logi FALSE
 - attr(*, "validate")= logi TRUE

Precision and Recall for FPReg

FPReg_df = provo_df %>% filter(metric == "FPReg") %>% spread(measure, value)
confusion_matrix <- table(FPReg_df$motr_value > 0, FPReg_df$eyetr_value > 0)
confusion_matrix

       
        FALSE TRUE
  FALSE    79 1518
  TRUE     22  918

true_positives <- confusion_matrix[2, 2]
false_positives <- confusion_matrix[2, 1]
false_negatives <- confusion_matrix[1, 2]

precision <- true_positives / (true_positives + false_positives)
recall <- true_positives / (true_positives + false_negatives)

print("precision of Motr FPReg:")

[1] "precision of Motr FPReg:"

print(precision)

[1] 0.9766

print("Recall of Motr FPReg:")

[1] "Recall of Motr FPReg:"

print(recall)

[1] 0.3768